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ABSTRACT 

Most research on the simulation of deformation and failure of metals has been and continues to be 
performed using the finite element method. However, the issues of mesh entanglement under large de- 
formation, considerable complexity in handling contact, and difficulties encountered while solving large 
deformation fluid-structure interaction problems have led to the exploration of alternative approaches. 
The material point method uses Lagrangian solid particles embedded in an Eulerian grid. Particles in- 
teract via the grid with other particles in the same body, with other solid bodies, and with fluids. Thus, 
the three issues mentioned in the context of finite element analysis are circumvented. 

In this paper, we present simulations of cylinders which fragment due to explosively expanding 
gases generated by reactions in a high energy material contained inside. The material point method is 
the numerical method chosen for these simulations discussed in this paper. The plastic deformation of 
metals is simulated using a hypoelastic-plastic stress update with radial return that assumes an additive 
decomposition of the rate of deformation tensor. Various plastic strain, plastic strain rate, and temper- 
ature dependent flow rules and yield conditions are investigated. Failure at individual material points 
is determined using porosity, damage and bifurcation conditions. Our models are validated using data 
from high strain rate impact experiments. It is concluded that the material point method possesses great 
potential for simulating high strain-rate, large deformation fluid-structure interaction problems. 

Keywords: Material Point Method, Fragmentation. 
INTRODUCTION 

The goal of this work is to present results from the simulation of the deformation and failure 
of a steel container that expands under the effect of gases produced by an explosively reacting 
high energy material (PBX 9501) contained inside. 

The high energy material reacts at temperatures of 450 K and above, This elevated tem- 
perature is achieved through external heating of the steel container. Experiments conducted at 
the University of Utah have shown that failure of the container can be due to ductile fracture 
associated with void coalescence and adiabatic shear bands. If shear bands dominate the steel 
container fragments, otherwise a few large cracks propagate along the cylinder and pop it open. 

Figure [T] shows the recovered parts of AISI 1026 steel containers after two different tests. 
The containers were initially 10 cm. in diameter and 0.6 cm thick. In the first test, shown in 
Figure [TJa), the container was heated over an open pool fire and shows ductile failure. In the 
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(a) Ductile fracture/Void Growth and Coalescence 





(b) Fragmentation/ Adiabatic Shear Bands 
FIG. 1. Experimental tests of exploding cylinders. 



second test, shown in Figure [jjb), the container was heated by means of electrical tape and 
fragmented after the explosion. 

The dynamics of the solid materials - steel and PBX 9501 - is modeled using the Lagrangian 
Material Point Method (MPM) ( |Sulsky et al. 1994[ ). Gases are generated from solid PBX 9501 



using a burn model (Long and Wight 2002). Gas-solid interaction is accomplished using an 
Implicit Continuous Eulerian (ICE) multi-material hydrodynamic code (Guilke y"et al. 2004] ). 
A single computational grid is used for all the materials. 



The constitutive response of PBX 9501 is modeled using ViscoSCRAM (Bennett et al. 
1998[ ), which is a five element generalized Maxwell model for the viscoelastic response cou- 
pled with statistical crack mechanics. Solid PBX 9501 is progressively converted into a gas 
with an appropriate equation of state. The temperature and pressure in the gas increase rapidly 
as the reaction continues. As a result, the steel container is pressurized, undergoes plastic 
deformation, and finally fragments. The entire process was simulated using the massively par- 



allel, Common Component Architecture (Armstrong et al. 1999 1 based, Uintah Computational 
Framework (UCF) ( |de St. Germain et al. 2000"1 ). 

The main issues regarding the constitutive modeling of the steel container are the selec- 
tion of appropriate models for nonlinear elasticity, plasticity, damage, loss of material stability, 
and failure. The numerical simulation of the steel container involves the choice of appropri- 
ate algorithms for the integration of balance laws and constitutive equations, as well as the 
methodology for fracture simulation. Models and simulation methods for the steel container 
are required to be temperature sensitive and valid for large distortions, large rotations, and a 
range of strain rates (quasistatic at the beginning of the simulation to approximately 10 6 s _1 at 
fracture). 

The approach chosen for the present work is to use hypoelastic-plastic constitutive models 
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that assume an additive decomposition of the rate of deformation tensor into elastic and plastic 
parts. Hypoelastic materials are known not to conserve energy in a loading-unloading cycle 
unless a very small time step is used. However, the choice of this model is justified under the 
assumption that elastic strains are expected to be small for the problem under consideration 
and unlikely to affect the computation significantly. 

Two plasticity models for flow stress are considered along with a two different yield con- 
ditions. Explicit fracture simulation is computationally expensive and prohibitive in the large 
simulations under consideration. The choice, therefore, has been to use damage models and 
stability criteria for the prediction of failure (at material points) and particle erosion for the 
simulation of fracture propagation. 

The outline of the paper is as follows. A brief description of the Material Point Method is 
given in Section [2j The stress update algorithm and the how various plasticity models, yield 
conditions, equations of state etc. are used during the stress update are discussed in Section [3] 
The models used for the simulations are discussed in Section[4] The results of some simulations 
are presented in Section [5] and conclusions are presented in Section [6] 

THE MATERIAL POINT METHOD 

The Material Point Method (MPM) ( |Sulsky et al. 1994[ ) is a particle method for structural 



mechanics simulations. In this method, the state variables of the material are described on 
Lagrangian particles or "material points". In addition, a regular, structured Eulerian grid is 
used as a computational scratch pad to compute spatial gradients and to solve the governing 
conservation equations. An explicit time-stepping version of the Material Point Method has 
been used in the simulations presented in this paper. The MPM algorithm is summarized 
below ( |Sulsky et al. 1995| ). 



It is assumed that an particle state at the beginning of a time step is known. The mass 
(to), external force (f ext ), and velocity (v) of the particles are interpolated to the grid using the 
relations 



id „ - > \ S gp nip , v 9 = ^2 S gp m p v p > f g Xt = ^2 S 9P ^ ^ 



p 



where the subscript (g) indicates a quantity at a grid node and a subscript (p) indicates a quantity 
on a particle. The symbol J2 P indicates a summation over all particles. The quantity (S gp ) is 
the interpolation function of node (g) evaluated at the position of particle (p). Details of the 
interpolants used can be found elsewhere (Bardenhagen and Kober 2004). 



Next, the velocity gradient at each particle is computed using the grid velocities using the 
relation 

VV P = Yl G 9P W 9 ( 2 ) 

g 

where G gp is the gradient of the shape function of node (g) evaluated at the position of particle 
(p). The velocity gradient at each particle is used to determine the Cauchy stress (<r p ) at the 
particle using a stress update algorithm. 

The internal force at the grid nodes (f^ nt ) is calculated from the divergence of the stress 
using 

ff = E G ^^ ( 3 ) 

p 

where V p is the particle volume. 
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The equation for the conservation of linear momentum is next solved on the grid. This 
equation can be cast in the form 

ui g a g = ^ l -f g nt (4) 

where a 9 is the acceleration vector at grid node (g). 

The velocity vector at node (g) is updated using an explicit (forward Euler) time integration, 
and the particle velocity and position are then updated using grid quantities. The relevant 
equations are 

V g (t+At)=Vg{t)+SLg At (5) 

v p (t + At) = v p (t) + ^2S gp a g At; x p (t + At) = x p (t) + ^ S gp v g At (6) 

a a 

The above sequence of steps is repeated for each time step. The above algorithm leads to 
particularly simple mechanisms for handling contact. Details of these contact algorithms can 



be found elsewhere (Bardenhagen et al. 2001 1. 



PLASTICITY AND FAILURE SIMULATION 

A hypoelastic-plastic, semi-implicit approach ( |Zocher et al. 2000 ) has been used for the 



stress update in the simulations presented in this paper. An additive decomposition of the rate 
of deformation tensor into elastic and plastic parts has been assumed. One advantage of this 
approach is that it can be used for both low and high strain rates. Another advantage is that 
many strain-rate and temperature-dependent plasticity and damage models are based on the 
assumption of additive decomposition of strain rates, making their implementation straightfor- 
ward. 

The stress update is performed in a co-rotational frame which is equivalent to using the 
Green-Naghdi objective stress rate. An incremental update of the rotation tensor is used instead 
of a direct polar decomposition of the deformation gradient. The accuracy of model is good if 
elastic strains are small compared to plastic strains and the material is not unloaded. It is also 
assumed that the stress tensor can be divided into a volumetric and a deviatoric component. The 
plasticity model is used to update only the deviatoric component of stress assuming isochoric 
behavior. The hydrostatic component of stress is updated using a solid equation of state. 

Since the material in the container may unload locally after fracture, the hypoelastic-plastic 
stress update may not work accurately under certain circumstances. An improvement would 
be to use a hyperelastic -plastic stress update algorithm. Also, the plasticity models are tem- 
perature dependent. Hence there is the issue of severe mesh dependence due to change of the 
governing equations from hyperbolic to elliptic in the softening regime ( |Hill and Hutchinson] 
1975} |Bazant and Belytschko 1985] Tvergaard and Needlema n 1990| ). Viscoplastic stress up- 



date models or nonlocal/gradient plasticity models (Ramaswamy and Aravas 1998a, Hao et al. 



2000) can be used to eliminate some of these effects and are currently under investigation. 

A particle is tagged as "failed" when its temperature is greater than the melting point of 
the material at the applied pressure. An additional condition for failure is when the porosity of 
a particle increases beyond a critical limit. A final condition for failure is when a bifurcation 
condition such as the Drucker stability postulate is satisfied. Upon failure, a particle is either 
removed from the computation by setting the stress to zero or is converted into a material with a 
different velocity field which interacts with the remaining particles via contact. Either approach 
leads to the simulation of a newly created surface. 
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Data: Persistent:Initial moduli, temperature, porosity, scalar damage, equation of state, 




plasticity model, yield condition, stability criterion, damage model 


Temporary: Particle state at time t 


Result: Particle state at time t + At 


for all the patches in the domain do 




Read the particle data and initialize updated data storage; 




for all the particles in the patch do 






Compute the velocity gradient, the rate of deformation tensor and the spin 






tensor; 






Compute the updated left stretch tensor, rotation tensor, and deformation 






gradient; 






Rotate the input Cauchy stress and the rate of deformation tensor to the material 






configuration; 






Compute the current shear modulus and melting temperature; 






Compute the pressure using the equation of state, update the hydrostatic stress, 






and compute the trial deviatoric stress; 






Compute the flow stress using the plasticity model; 






Evaluate the yield function; 






if particle is elastic then 








Rotate the stress back to laboratory coordinates; 








Update the particle state; 






else 








Find derivatives of the yield function; 








Do radial return adjustment of deviatoric stress; 








Compute updated porosity, scalar damage, and temperature increase due to 








plastic work; 








Compute elastic-plastic tangent modulus and evaluate stability condition; 








Rotate the stress back to laboratory coordinates; 








Update the particle state; 








if Temperature > Melt Temperature or Porosity > Critical Porosity or 








Unstable then 








Tag particle as failed; 








end 






end 




end 




end 






Convert failed particles into a material with a different velocity field; 



Algorithm 1: Stress Update Algorithm 



In the parallel implementation of the stress update algorithm, sockets have been added 
to allow for the incorporation of a variety of plasticity, damage, yield, and bifurcation models 
without requiring any change in the stress update code. The algorithm is shown in Algorithm[T] 
The equation of state, plasticity model, yield condition, damage model, and the stability crite- 
rion are all polymorphic objects created using a factory idiom in C++ (Coplien 1992 ). 

MODELS 
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The stress in the solid is partitioned into a volumetric part and a deviatoric part. Only 
the deviatoric part of stress is used in the plasticity calculations assuming isoschoric plastic 
behavior. 

The hydrostatic pressure (p) is calculated either using the bulk modulus (K) and shear mod- 



ulus (p) or from a temperature-corrected Mie-Gruneisen equation of state of the form (Zocher 



etal. 2000) 



P 



Po cic [i + (i - 30 c] 



fl 



+ T C P T 



( = (p/po - 1) 



(7) 



(Sa ~ 1)C] 

where Co is the bulk speed of sound, po is the initial density, p is the current density, C p is 
the specific heat at constant volume, T is the temperature, To is the Gruneisen's gamma at 
reference state, and S a is the linear Hugoniot slope coefficient. 

Depending on the plasticity model being used, the pressure and temperature-dependent 
shear modulus (p) and the pressure-dependent melt temperature (T m ) are calculated using the 



relations (Steinberg et al. 1980) 



P = Mo 



1 + A 



P 

,1/3 



B(T - 300) 



l m0 exp 



2a 1 



2(T -o-l/3) 



(8) 
(9) 



where, po is the shear modulus at the reference state(T = 300 K, p = 0, e p = 0), e p is the plastic 
strain, rj = p/po is the compression, A = (1/ po)(dp/dp), B = (1/ po)(dp/dT), T m o is the 
melt temperature at p = po, and a is the coefficient of the first order volume correction to 
Gruneisen's gamma. 

We have explored two temperature and strain rate dependent plasticity models - the Johnson- 



Cook plasticity model (Johnson and Cook 1983) and the Mechanical Threshold Stress (MTS) 



plasticity model (Follansbee and Kocks 1988 ; Goto et al. 2000 1. The flow stress (oy) from the 
Johnson-Cook model is given by 



a f = [A + B(e P T][l + Cln(e*)][l - (T*) m ] ; e* 



£p0 



(T 



(10) 



where e p o is a user defined plastic strain rate, A, B, C, n, m are material constants, T r is the 
room temperature, and T m is the melt temperature. 
The flow stress for the MTS model is given by 



where 



M = Mo 



Si 



D 

exp (f) - 1 
kT e 0i 



a a + —Si&i + —S e a e 
Po Po 



(11) 
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e = e [l-F(X)] + e IV F(X) ; 9 = a + ai In e + a 2 Vl-a 3 T 



X 



Pes 

& (n+l) = a (n) + 0Ae 



; F(X) = tanh(aX) ; ln(<7 es /<7, 



esOj 



\nb 3 g< 



In - — 
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and cr a is the athermal component of mechanical threshold stress, /iq is the shear modulus at 
K, D, To are empirical constants, &i represents the stress due to intrinsic barriers to thermally 
activated dislocation motion and dislocation-dislocation interactions, <r e represents the stress 
due to microstructural evolution with increasing deformation, k is the Boltzmann constant, b is 
the length of the Burger's vector, <7o[i,e] w& tne normalized activation energies, £o[j,e] are con- 
stant strain rates, qu e }iP\i,e] we constants, 9q is the hardening due to dislocation accumulation, 
ao, ai, a 2 , 03, Oiv, ct are constants, a es is the stress at zero strain hardening rate, a es o is the 
saturation threshold stress for deformation at K, goes is a constant, and e es o is the maximum 
strain rate. 

We have decided to focus on ductile failure of the steel container. Accordingly, two yield 
criteria have been explored - the von Mises condition and the Gurson-Tvergaard-Needleman 
(GTN) yield condition ( Gurson 1977[ Tvergaard and Needleman 1984 1 which depends on 
porosity. An associated flow rule is used to determine the plastic rate parameter in either case. 
The von Mises yield condition is given by 



"eg 



1 = 0; a eq 



(12) 



where a eq is the von Mises equivalent stress, a d is the deviatoric part of the Cauchy stress, and 
is the flow stress. The GTN yield condition can be written as 



- ei ) + 2qi /* cosh ( q 2 — 



Tr(a) 



(1 + qzft) = 



(13) 



where qi,q 2 , q% are material constants and is the porosity (damage) function given by 



fc + Kf-fc 



for / < f c 
for f > f c 



(14) 



where is a constant and / is the porosity (void volume fraction). The flow stress in the matrix 
material is computed using either of the two plasticity models discussed earlier. Note that the 
flow stress in the matrix material also remains on the undamaged matrix yield surface and uses 
an associated flow rule. 

The evolution of porosity is calculated as the sum of the rate of growth and the rate of 



nucleation (Ramaswamy and Aravas 1998b I. The rate of growth of porosity and the void 



nucleation rate are given by the following equations (Chu and Needleman 1980) 



/ — /nucl ~\~ /grow 

/grow = (1 - /)Tr(D p ) 



/, 



fn 



nucl 



exp 



1 (e p - e n ) s 

2 & 



(15) 
(16) 

(17) 
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where D p is the rate of plastic deformation tensor, /„ is the volume fraction of void nucleat- 
ing particles , e n is the mean of the distribution of nucleation strains, and s n is the standard 
deviation of the distribution. 

Part of the plastic work done is converted into heat and used to update the temperature of a 
particle. The increase in temperature (AT) due to an increment in plastic strain (Ae p ) is given 



by the equation (Borvik et al. 2001 ) 



AT 



P C P P 



(18) 



where x is the Taylor-Quinney coefficient, and C p is the specific heat. A special equation for 



the dependence of C p upon temperature is also used for steel (Goto et al. 2000). 



C p = 10 3 (0. 09278 + 7.454 x 10~ 4 T + 12404.0/T 2 



(19) 



Under normal conditions, the heat generated at a material point is conducted away at the 
end of a time step using the heat equation. If special adiabatic conditions apply (such as in 
impact problems), the heat is accumulated at a material point and is not conducted to the 
surrounding particles. This localized heating can be used to simulate adiabatic shear band 
formation. 

After the stress state has been determined on the basis of the yield condition and the as- 
sociated flow rule, a scalar damage state in each material point can be calculated using either 



of two damage models - the Johnson-Cook model (Johnson and Cook 1985 ) or the Hancock- 



MacKenzie model (Hancock and MacKenzie 1976). While the Johnson-Cook model has an 



explicit dependence on temperature, the Hancock-McKenzie model depends on the tempera- 
ture implicitly, via the stress state. Both models depend on the strain rate to determine the value 
of the scalar damage parameter. 

The damage evolution rule for the Johnson-Cook damage model can be written as 



D 



J 



D\ + T>2 exp 



^3 



-a 



[1 + Z) 4 ln(ep*)][l + D 5 r*] ; a* 



Tr(<x) 



eq 



(20) 

where D is the damage variable which has a value of for virgin material and a value of 1 at 
fracture, e p is the fracture strain, D\, D2, -D3, TJ4, D5 are constants, a is the Cauchy stress, and 
T* is the scaled temperature as in the Johnson-Cook plasticity model. 
The Hancock-Mac Kenzie damage evolution rule can be written as 



D 



J 



1.65 



exp(1.5<7*) 



(21) 



The determination of whether a particle has failed can be made on the basis of either or all 
of the following conditions: 

• The particle temperature exceeds the melting temperature. 



The TEPLA-F fracture condition (Johnson and Addessio 1988) is satisfied. This condi- 
tion can be written as 

(f/fcf + (e P /elf = 1 (22) 

where / is the current porosity, f c is the maximum allowable porosity, e p is the current 
plastic strain, and e p is the plastic strain at fracture. 
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An alternative to ad-hoc damage criteria is to use the concept of bifurcation to deter- 
mine whether a particle has failed or not. Two stability criteria have been explored in 
this paper - the Drucker stability postulate ( jDruck er 1959 ) and the loss of hyperbolic- 
ity criterion (using the determinant of the acoustic tensor) ( [Rudnicki and Rice 1975[ 



Perzyna 1998). 



The simplest criterion that can be used is the Drucker stability postulate ( |Drucker 1959 1 



which states that time rate of change of the rate of work done by a material cannot be negative. 
Therefore, the material is assumed to become unstable (and a particle fails) when 

& : D p < (23) 

Another stability criterion that is less restrictive is the acoustic tensor criterion which states 
that the material loses stability if the determinant of the acoustic tensor changes sign (Rud 



nicki and Rice 1975] Perzyna 1998 1. Determination of the acoustic tensor requires a search 
for a normal vector around the material point and is therefore computationally expensive. A 
simplification of this criterion is a check which assumes that the direction of instability lies in 



the plane of the maximum and minimum principal stress (Becker 2002 1. In this approach, we 
assume that the strain is localized in a band with normal n, and the magnitude of the velocity 
difference across the band is g. Then the bifurcation condition leads to the relation 

R ij3j = i R ij = M ikjinkni + M Ukj n k ni - a ik n^n k (24) 

where M; L j k i are the components of the co-rotational tangent modulus tensor and are the 
components of the co-rotational stress tensor. If det(Rij) < 0, then gj can be arbitrary and 
there is a possibility of strain localization. If this condition for loss of hyperbolicity is met, 
then a particle deforms in an unstable manner and failure can be assumed to have occurred at 
that particle. 

SIMULATIONS 

The first set of simulations was performed using the geometry shown in Figure^a). A steel 
cylinder was used to confine the PBX 9501 material and the simulation was started with both 
materials at a temperature of 600 K. At this temperature, PBX 9501 reacts and forms gases 
which expand the cylinder. A quarter of the cylinder was modeled using a 160 x 160 x 1 grid 
with 8 particles per grid cell. The shapes of the cylinder after failure for two different materials 
are shown in Figure [2] 

The simulation shown in Figure |2jb) was performed using material data for 4340 steel, a 
Mie-Griineisen equation of state, the Johnson-Cook flow stress model, the Gurson yield condi- 
tion, the Johnson-Cook damage model, and checks of both the Drucker stability postulate and 
the loss of hyperbolicity condition. The simulation of a HY 100 steel container shown in Fig- 
ure [2jc) was performed using a Mie-Griineisen equation of state, the MTS flow stress model, 
the Gurson yield condition, the Hancock-MacKenzie damage model, and the same stability 
checks as the 4340 steel. The material properties and the parameters used in the models are 
shown in Table [TJ The materials are given an initial mean porosity of 0.005 using a Gaussian 
distribution with a standard deviation of 0.001 and a mean scalar damage value of 0.01 with a 
standard deviation of 0.005. 

The expected number of fragments (N) along the circumference of the exploding cylinder 
can be approximated using the following analytical result (Gra dy and Hightower 1992) 



N = 2 ^ Pj ^f) 1/3 (25) 
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(a) Geometry 




(b) 4340 Steel. (c) HY 100 Steel. 

FIG. 2. Simulations of fragmenting cylinders - two-dimensional view. 

where p is the density, Ro is the initial cylinder radius, V is the expansion velocity at the radius 
of fracture, and F is the fragmentation energy. 

The fragmentation energy in tension (IV) and in shear (F$) are given by 

2 1, a \ Y 6 a ZA f I 

where K c is the fracture toughness, E is the Young's modulus, p is the density, C p is the specific 
heat at constant pressure, a is the thermal softening coefficient, x i s tne thermal diffusion 
coefficient, Y is the yield strength in simple tension, and 7 is the shear strain rate. 

For the expanding 4340 steel cylinder of that we have simulated, the relevant quantities 
are p = 7830 kg/m 3 , E = 208 GPa, K c = 80 MN/m 2 / 3 , Y = 792 MPa, C v = All J/kg K, x = 
1.5xl0~ 5 m 2 /s, a = 7.5xl0" 4 /K, R = 0.054 m, V = 300 m/s, 7 = 1000 /s, Ft = 1.5 xlO 4 
J/m 2 , and F$ = 5.2 xlO 4 J/m 2 . Accordingly, the expected number of fragments (for the whole 
cylinder) are iV (tension) = 29 and N (shear) = 20. For a quarter of the cylinder, the number 
of fragments is expected to be between 8 and 5. We get approximately 6 to 7 fragments in our 
simulations, which implies that our results are qualitatively acceptable. Both the steels show 
similar fragmentation though the exact shape of the fragments differs slightly. For this reason, 
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TABLE 1. Material Properties and Parameters for Steels. 



4340 Steel properties and Johnson-Cook parameters 


P 




T 


K 


f 


X 










(kg/m 3 ) 


(MPa m 3 /kj 


; K) (K) 


(GPa) 


(GPa) 












7830.0 


477.0 


1793.0 


173.3 


80.0 


0.9 










A 


B 


C 


n 


m 


Dx 


D 2 


D 3 






(MPa) 


(MPa) 


















792.0 


510.0 


0.014 


0.26 


1.03 


0.05 


3.44 


-2.12 


0.002 


0.61 


HY100 Steel properties and MTS parameters 


P 


C p 


T 




n 


X 










(kg/m 3 ) 


(MPa m 3 /kf 


! K) (K) 


(GPa) 


(GPa) 












7860.0 


477.0 


2000.0 


150.0 


69.0 


0.9 












Mo 


D 


To 


k/b 3 


m 










(MPa) 


(GPa) 


(GPa) 


(K) 


(xlO 6 ) 






(xlO 13 ) 


(xlO 7 ) 




40.0 


71.46 


2.9 


204 


0.905 


1.161 


1.6 


1.0 


1.0 




Pi 


Qi 


Pe 


9e 


&i 
(MPa) 


ao 
(xlO 9 ) 


a i 


a 2 


(xlO 6 ) 




0.5 


1.5 


0.67 


1.0 


1341 


6 








2.0758 




9 iv 


a 


CesO 




O'esO 












(xlO 6 ) 




(xlO 7 ) 




(MPa) 












200.0 


3 


1.0 


0.112 


822.0 












Mie-Gruneisen equation of state parameters 


Co 


To 


See 
















(m/s) 




















3574 


1.69 


1.92 
















GTN yield condition and porosity evolution parameters 




12 






fc 


Jn 


s« 








1.5 


1.0 


2.25 


4.0 


0.05 


0.1 


0.3 


0.1 







the three-dimensional simulations were performed using 4340 steel and the associated models 
discussed above. 

Figure [3] shows the fragmentation obtained from three-dimensional simulations of a cylin- 
der with end-caps. A quarter of the geometry is modeled, assuming symmetry. The cylinder 
is made of 4340 steel and contains PBX 9501. The simulation is started with both materials 
at a temperature of 600 K. A hypoelastic constitutive model is used to determine the volumet- 
ric response of the material. The Johnson-Cook plasticity model is used to calculate the flow 
stress. The von Mises yield condition is used to determine the boundary of the elastic and 
plastic domains. A Johnson-Cook damage model is used to compute a scalar damage param- 
eter. A uniform initial porosity is assigned to all steel particles and evolved according to the 
models discussed in the previous section. A particle is deemed to have failed when the modi- 
fied TEPLA-F condition is satisfied, the temperature is more than the melting temperature, or 
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(a) Fragments of the container. (b) Gases escaping from the container. 

FIG. 3. Simulations of fragmenting cylinders - three-dimensional view. 

the Drucker stability postulate/loss of hyperbolicity condition is satisfied. Upon failure, the 
particle stress is set to zero. 

The simulations capture some of the qualitative features observed in the experiments of 
steel cylinders heated using heat tapes. Some high particle velocities are observed upon failure. 
Simulations have shown that these velocities are due to some increase in the total energy due to 
the setting of the particle stress to zero. Computations where failed particles are converted into 
a material with a different velocity field are currently under way along with other validation 
efforts to quantify the error in the calculations. 

Simulations have also been performed on containers heated by a pool fire. Four snapshots 
of one such simulation for 4340 steel using the Johnson-Cook plasticity and damage models, a 
Mie-Gruneisen equation of state, the von Mises yield condition and uniform initial porosity are 
shown in Figure [4] The fire is simulated by a hot jet of air and interacts with the container via 
heat conduction and momentum transport. After some time, the contents of the container reach 
ignition temperature and reaction proceeds rapidly. Axial cracks open up in the container that 
are qualitatively similar to those observed in experiments. These cracks join to form fragments 
which interact with the fire. More details of such simulations can be found elsewhere (Guilkey 
et al. 2004| ). Validation of the pool fire and container interaction scenario is currently being 
performed in collaboration with researchers at the Lawrence Livermore National Laboratory. 

SUMMARY AND CONCLUSIONS 

A computational scheme for the simulation of the fragmentation of cylinders due to inter- 
action with gases from a reacting high energy material has been presented. The scheme allows 
for the incorporation of various plasticity models, yield conditions, damage models, equations 
of state, and stability checks within the same stress update code. A number of such models 
have been listed and the corresponding material properties and parameters for two steels have 
been collected from various sources and presented in a compact form. 

Simulations of exploding cylinders in two-dimensions have been compared with analytical 
solutions for the expected number of fragments and found to provide qualitative agreement. 
Three-dimensional simulations also show qualitative agreement with experiments in the direc- 
tions of the dominant cracks. Snapshots from the simulation of a fully coupled container-fire 
simulation have also shown qualitative agreement with experiments. 
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FIG. 4. Simulations of fragmentation of a cylinder heated by a fire. 

Two issues that have been identified as important for the simulation are the conservation of 
energy and mesh dependence of the results. Validation simulations that are currently underway 
have shown that energy is better conserved when particles are converted into a material with a 
different velocity field after failure (rather than setting the stress to zero upon failure). Results 
of these tests will be presented in future work. In the absence of a limiting length scale in the 
computation, strongly mesh dependent behavior can be expected in the softening regime of the 
stress-strain relationship. This mesh dependence occurs when we use temperature dependent 
elastic/plastic constitutive equations and when we degrade the yield strength of the porous 
material using porosity. One way of minimizing mesh dependence is to use a rate-dependent 
stress update algorithm. However, such an approach is not sufficient to remove the effects of 
all the possible causes of mesh dependence. Validation experiments are currently underway to 
determine the extent of mesh dependence and nonlocal/gradient plasticity approaches that can 
be formulated for the material point method. Overall, the material point method appears to be a 
promising approach for simulating high rate, coupled fluid- structure interaction problems and 
fragmentation. 
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